Create figures for the clinical paper

All data is coming from tables for the paper

# labels for time points
timeLab  = c("Baseline", "Post Run-in", "Week 8")

# set colors
recistCol = c("Non-responder" = "red","Responder" = "blue")
cbrCol = c("CR,PR,SD" = "blue","PD" = "red")
subtypeCol = c("HR+" = "red", "TNBC" = "blue")
timeCol = setNames(c("black", 'grey45', 'grey75'), timeLab)

# set root folder
rootDir = "./"

# Load data patient data
dta_pts = read.csv(paste0(rootDir,"input_tables/table_24patient_info.csv"), row.names = 1)
dta_pts$ID = rownames(dta_pts)

# remove patients with no response
dta_pts = dta_pts %>% filter(Best_Resp != "Unknown")


#==========================
# file with TILs and PDL1 data
cytof_dat = read.csv(paste0(rootDir,"input_tables/table_cytof_marker_expression.csv"))

#markers = c("CD40", "CD103", "CCR5", "CD141")
markers = colnames(cytof_dat)[4:ncol(cytof_dat)]
#===========================
# add patient info 
cytof_dat = cbind(cytof_dat,  dta_pts[cytof_dat$patientID,])

# take only baseline and time 1 values
df_BvT1 <- cytof_dat %>% filter(timePoint %in% c('Baseline','Post Run-in'))

# take only time 1 and time 2 values
df_T1vT2 <- cytof_dat %>% filter(timePoint %in% c('Week 8','Post Run-in'))

P-values

Pairwise Wilcoxon test was used for each time point comparing responders vs non-responders. And paired Wilcoxon test was used for baseline vs post run-in and post run-in vs week 8. We took the differences between all time points and then compare those differences between responders vs non-responders to see if there is significant difference in changes after treatment

# Calculate p-values from Wilcoxon test: paired for time points and not paired for response
res = c()
# split by clusters
cytof_dat_byClust = split(cytof_dat, f = factor(cytof_dat$cluster))
for(j in names(cytof_dat_byClust))
{
  #
  d = cytof_dat_byClust[[j]]
  for (k in markers)
  {
    # create per patient matrix with three time points in columns
    # some values would be NA
    mat <- data.frame(matrix(nrow = length(unique(d$patientID)),
      ncol = 3, dimnames = list(unique(d$patientID),
                      levels(factor(d$timePoint)))), check.names = F)
    # fill in matrix
    for(i in levels(factor(d$timePoint)))
    {
      d1 = d %>% filter(timePoint == i)
      mat[d1$patientID, i] = d1[,k]
    }
    # add response
    mat$RECIST = cytof_dat[rownames(mat),"RECIST"]
    # create treatment effect 
    # substract baseline from post run-in and week 8 
    mat$delta_T1_B =  mat[,2]-mat[,1]
    mat$delta_T2_B =  mat[,3]-mat[,1]
    # and week 8 from post run-in 
    mat$delta_T2_T1 =  mat[,3]-mat[,2]
    
    # running paired test
    # Baseline vs post run-in
    p1 = wilcox.test(mat[,1], mat[,2], paired = T)$p.value
    # post run-in vs week 8
    p2 = wilcox.test(mat[,2], mat[,3], paired = T)$p.value
    
    # test by response
    # for baseline 
    p3 = wilcox.test(Baseline ~ RECIST, data = mat)$p.value
    # for post run-in 
    p4 = wilcox.test(mat[,2] ~ mat[,"RECIST"], data = mat)$p.value
    # for week 8
    p5 = wilcox.test(mat[,3] ~ mat[,"RECIST"], data = mat)$p.value
    # compare differences in treatment effect in responders vs non-responders
    p6 = wilcox.test(delta_T1_B ~ RECIST, data = mat)$p.value
    p7 = wilcox.test(delta_T2_B ~ RECIST, data = mat)$p.value
    p8 = wilcox.test(delta_T2_T1 ~ RECIST, data = mat)$p.value
   
    # combine all p-values and add means
    res = rbind(res, c(marker = k, cluster = j, BvsT1 = round(p1,3),
                       T1vsT2 = round(p2,3), B_RvsNR = round(p3,3), 
                       T1_RvsNR = round(p4,3), T2_RvsNR = round(p5,3), 
                       d_T1_B_RvsNR = round(p6,3), 
                       d_T2_B_RvsNR = round(p7,3), 
                       d_T2_T1_RvsNR = round(p8,3),
                       mean_d_T1_B_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T1_B) %>% unlist(), na.rm = T),3),
                       mean_d_T1_B_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T1_B) %>% unlist(), na.rm = T),3),
                       mean_d_T2_B_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T2_B) %>% unlist(), na.rm = T),3),
                       mean_d_T2_B_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T2_B) %>% unlist(), na.rm = T),3),
                       mean_d_T2_T1_R = round(mean(mat %>% filter(RECIST == "Responder") %>% select(delta_T2_T1) %>% unlist(), na.rm = T),3),
                       mean_d_T2_T1_NR = round(mean(mat %>% filter(RECIST == "Non-responder") %>% select(delta_T2_T1) %>% unlist(), na.rm = T),3)))
  }
}
res = res[order(res[,'marker']),]

datatable(res)
write.csv(res,file = "table_cytof_wilcox_all_markers.csv", row.names = F)


# make boxplots for markers/clusters that have p-value < 0.05 in d_T1_B_RvsNR
sets = data.frame(res)
sets = sets %>% filter(d_T1_B_RvsNR < 0.05 & cluster%in% c('B','CD4_T','CD8_T','DNT','DPT','NK_like',"cDC1"))

for (i in 1:nrow(sets))
{
  k = sets[i,'cluster']
  j = sets[i,'marker']
d = cytof_dat %>% filter(cluster == k)
  df = cbind(protein = d[,j], 
             d[,c("patientID", "timePoint", "RECIST","Subtype") ] )

    p = ggpaired(df , x="timePoint",y="protein", id = "patientID", 
                 line.size = 0.4, title = paste0(j,". ", k), 
                 line.color = "RECIST",palette =  recistCol, xlab="Time") +
      facet_wrap("RECIST") +
      scale_x_discrete(labels=timeLab)
    
    print(p)
}

# additionally, plot CD103 in cDC1 cluster

  k = "cDC1"
  j = "CD103"
d = cytof_dat %>% filter(cluster == k)
d = CreateAllFacet(d,"RECIST")
  df = cbind(protein = d[,j], 
             d[,c("patientID", "timePoint", "RECIST","Subtype", "facet") ] )

    p = ggpaired(df , x="timePoint",y="protein", id = "patientID", 
                 line.size = 0.4, title = paste0(j,". ", k), 
                 line.color = "RECIST",palette =  recistCol, xlab="Time") +
      facet_wrap("facet") +
      scale_x_discrete(labels=timeLab)
    
    print(p)

Baseline vs Post Run-in

#ggpaired(df_BvT1 %>% filter(cluster %in% c("GMDSC")),
for (i in markers)
{
  p = ggpaired(df_BvT1, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
   palette =  recistCol, title = paste0(i,". Baseline vs Post Run-in"), xlab="Time", ylab="MMI") + 
    scale_x_discrete(labels=timeLab)+ 
    facet_wrap("cluster", scales = 'free') +
    theme(axis.text.x=element_text(angle = 45, hjust = 1))

  print(p)
}

Post Run-in vs Week 8

for (i in markers)
{
  p = ggpaired(df_T1vT2, x="timePoint",y=i, id = "patientID", line.color = "RECIST", line.size = 1,
   palette =  recistCol, title = paste0(i,". Post Run-in vs Week 8"), xlab="Time", ylab="MMI") +
#    scale_x_discrete(labels=timeLab)+
    facet_wrap("cluster", scales = 'free') +
    theme(axis.text.x=element_text(angle = 45, hjust = 1))

  print(p)
}

sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] scales_1.2.1      gridExtra_2.3     DT_0.28           data.table_1.14.8
## [5] ggrepel_0.9.3     ggpubr_0.6.0      ggplot2_3.4.2     dplyr_1.1.2      
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.7        utf8_1.2.3        generics_0.1.3    tidyr_1.3.0      
##  [5] rstatix_0.7.2     digest_0.6.33     magrittr_2.0.3    evaluate_0.21    
##  [9] grid_4.3.1        fastmap_1.1.1     jsonlite_1.8.7    backports_1.4.1  
## [13] purrr_1.0.1       fansi_1.0.4       crosstalk_1.2.0   jquerylib_0.1.4  
## [17] abind_1.4-5       cli_3.6.1         rlang_1.1.1       ellipsis_0.3.2   
## [21] munsell_0.5.0     withr_2.5.0       cachem_1.0.8      yaml_2.3.7       
## [25] tools_4.3.1       ggsignif_0.6.4    colorspace_2.1-0  broom_1.0.5      
## [29] vctrs_0.6.3       R6_2.5.1          lifecycle_1.0.3   car_3.1-2        
## [33] htmlwidgets_1.6.2 pkgconfig_2.0.3   pillar_1.9.0      bslib_0.5.0      
## [37] gtable_0.3.3      glue_1.6.2        Rcpp_1.0.11       highr_0.10       
## [41] xfun_0.39         tibble_3.2.1      tidyselect_1.2.0  rstudioapi_0.15.0
## [45] knitr_1.43        farver_2.1.1      htmltools_0.5.5   labeling_0.4.2   
## [49] rmarkdown_2.23    carData_3.0-5     compiler_4.3.1